Mon. Not. R. Astron. Soc. 000, [T]-?? (2008) Printed 17 November 2011 (MN lAT^ style file v2.2) 



Global structure and kinematics of stellar haloes in 
cosmological hydrodynamic simulations 

I. G. McCarthyi'2*, A. S. Font^, R. A. Crain^'^, A. J. Deason^, J. Schaye^ 
T. Theuns^'^ 

^ Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CBS OH A 
'^Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CBS OHA 

Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, Victoria SI 22, Australia 
'^Leiden Observatory, Leiden University, P. O. Box 951 S, 2300 RA Leiden, the Netherlands 

^Institute of Computational Cosmology, Department of Physics, University of Durham, Science Laboratories, South Road, Durham DHl SLE 
^Department of Physics, University of Antwerp, Campus Groenenborger, Groenenborgerlaan 171, B-2020 Antwerp, Belgium 

Accepted ... Received ... 

ABSTRACT 

Wc use the Galaxies-Intergalactic Medium Interaction Calculation (GIMIC) suite of 
cosmological hydrodynamical simulations to study the global structure and kinemat- 
ics of stellar spheroids of Milky Way mass disc galaxies. Font et al. have recently 
demonstrated that these simulations are able to successfully reproduce the satellite 
luminosity functions and the metallicity and surface brightness profiles of the spheroids 
of the Milky Way and M31. A key to the success of the simulations is a significant 
contribution to the spheroid from stars that formed in situ. While the outer halo is 
dominated by accreted stars, stars formed in the main progenitor of the galaxy dom- 
inate at r < 30 kpc. In the present study we show that this component was primarily 
formed in a proto-disc at high redshift and was subsequently liberated from the disc 
by dynamical heating associated with mass accretion. As a consequence of its origin, 
the in situ component of the spheroid has different kinematics (namely net prograde 
rotation with respect to the disc) than that of the spheroid component built from the 
disruption of satellites. In addition, the in situ component has a flattened distribution, 
that is due in part to its rotation. We make comparisons with measurements of the 
shape and kinematics of local galaxies, including the Milky Way and M31, and stacked 
observations of more distant galaxies. We find that the simulated disc galaxies have 
spheroids of the correct shape (oblate with a median axis ratio of ~ 0.6 at radii of 
< 30 kpc, but note there is significant system-to-system scatter in this quantity) and 
that the kinematics show evidence for two components (due to in situ vs. accreted), 
as observed. Our findings therefore add considerable weight to the importance of dis- 
sipative processes in the formation of stellar haloes and to the notion of a 'dual stellar 
halo'. 

Key virords: Galaxy: evolution — Galaxy: formation — Galaxy: halo — galaxies: 
evolution — galaxies: formation — galaxies: haloes 



1 INTRODUCTION 

It has long been recognised that the extended, diffuse stellar 
haloes of galaxies contain a potential treasure trove of infor- 
mation about the formation and evolution of galaxies (e.g . , 
lEggen. Lvnden-Bell fc Sandagd [l96^ : ISearle fc ZiiinlllQTSh . 
Observationally, we are entering a 'golden age' for stud- 
ies of stellar haloes, with the exploitation of large surveys 
such as the SDSS / Sloan Extension for Galactic Under- 

* E-mail: mccarthy@ast.cam.ac.uk 



standing and Exploratioi£] (SEGUE, lYannv et all 120091 ') . 
RAVeQ, SkyMappeiQ, Pan-Starrfl and Gai43that are yield- 
ing (or will yield) exquisite structural and kinematic mea- 
surements of the Milky Way's stellar halo and its sub- 
structures. At the same time, dedicated surveys of nearby 

http://www.sdss3.org/ 
^ http:/ /www. rave-survey.aip.de/rave/ 
^ http:/ /www. mso.anu.edu.au/skymapper/ 
* http:/ /pan-starrs. ifa.hawaii.edu/public/ 
^ http://gaia.esa.int/ 
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galcixies, such as M31 (e.g ., Guhathakurta et al. _ 20051: 



Irwin et alj|2005l:llbata et alj|2007l ). NGC 891 (jlbata et all 
20091 : iMouhcine et al.ll2oiol V and other nearby galaxies (e.g.. 



in the GHOSTS survejp) , provide us with a valuable exter- 
nal perspective and a point of comparison for the Milky Way 
(i.e., how typical is the Milky Way?). 

Over the past few years there have also been signifi- 
cant advances in the modelling of stellar haloes of galaxies 
with masses similar to that of the Milky Way. The major- 
ity of these studies have relied on so-called 'hybrid' models, 
which couple simple recipes for assigning stellar mass to in- 
falling satellites in dissipationless (i.e., dark matter only) 
cosmological simulat i ons (e .g.. Bullock & Johnst onI l2005l: 
Dicman d et~aLl l2005l: iFont et all l2006al : lOe Lucia fc Helmil 



2008. : .Cooper et al.ll201Cl ). In terms of assigning stellar mass 
to the infalling satellites, the models are generally tuned 
to reproduce the observed luminosity function of satellites 
around the Milky Way. The stellar mass is distributed by 
tagging centrally-located dark matter particles and track- 
ing the DM particles using their IDs. In general, the mod- 
els have been largely successful at reproducing the observed 
outer stellar halo, as well as the observed demographics of 
stellar streams. 

However, it has also become clear that such h ybrid mod- 



els hav e a number of problems. For example. Font et al.l 



2006bh . lDe Lucia fc Helmil ()2008h . and lCooper et al.1 ()201GI 1 



have all shown that they fail to produce significant metal- 
licity gradients in the stellar s pheroids, whereas o b serva- 
tions of the Milky Way (e. g.. ICaroUo et al.1 l200l I2OI0I : 



i Kalirai et al 



2006 



Ide Jong etU] |2010| . but see Sesar et al.' l201 l|) and M31 



iKoch et aLr 2008: Gilber t et al.1 12009| ) 



show evidence for large-scale gradie ntfl Furthermo r e, ob - 



servations of t he Milky Way (e.g.. IChiba fc BeersI I2OOII ') 



and M31 (e.g.. iPritchet fc van den Berg hj 1994), as well a s 



stacked observations of distant galaxies I zibetti et al.ll2004h . 
show that the spheroid has a highly flattened (oblate) dis- 



^ http://www.stsci.edu/ djrs/ghosts/ 

^ There is a degree of confusion in the literature over the term 
"metallicitv gradient". Detailed observations of main-sequence 



turn-off stars fe.a:.. ICarolIo et al.ll2007l. 


2OIOI: Ide Jo 


e et al.ll2010l 


and blue horizontal-branch stars (e.g.. 


Beers et al. 


2OIII) in the 



Milky Way show evidence for two components (one relatively 
metal rich and the other relatively metal poor) in the metal- 
licity distribution function (MDF) that overlap spatially in the 
inner regions of the halo. In the outer regions, beyond r ~ 20 
kpc or so, however, the metal-poor component is dominant and 
there is no further alteration of the MDF as one progresses further 
into the outer halo (i.e., no gradient). The transition from metal- 
rich component dominance to metal-poor component dominance 
would suggest a gradient if these populations are not separated 
explicitly (e.g., as in the case of observations of M31), because of 
the lessening degree of dominance between the metal-rich com- 
ponent, relative to the metal-poor component, as one moves out- 
ward. It is in this sense that we use the term gradient. iFont et al.l 
demonstrated that the simulated galaxies in GIMIC exhibit 
exactly this behavior; i.e., that there are two populations (in situ 
and accreted) with different characteristic metallicities that over- 
lap spatially, and it is the change in the mass fraction of these 
populations with radius that gives rise to an overall gradient out 
to radii of 30 kpc (see their Fig. 8). Beyond this radius, the ac- 
creted component dominates by mass, thus there is no further 
decline in the mean metallicity. 



tribution, whereas the stellar haloes in hybrid simulations 
usually have a pro late configuration like the dark matter 
ijCooper et al.ll2010l ). Observations of the spheroid in the So- 
lar neighbourhood show kinematic and chemodynamic evi- 
dence for two distinct components to the spheroid (the 'dua l 
halo' e.g.. lCarollo et al ]|200l l2010l : iNissen fc Schusteill2010l ) 
that may be difficult to reconcile with purely dissipationless 
hybrid methods. Finally, hybrid models appear to overpro- 
duce the degree of s ubstructure seen in the stellar hal o of 
the Milky Way (e.g.. lHelmi et al.ll201ll : IXue et al.ll201ll ). 

Hybrid models implicitly assume that the spheroid 
is built entirely of the tidal disruption of infalling satel- 
lites, they do not include/allow for dissipational processes 
such as in situ star formation (i.e., star formation in 
the spheroid itself, or star formation in a disc which 
is subsequently heated to become a spheroid). Studies 
of the formation of stellar haloes in cosmological hydro- 
dynamic simulations, however, routinely show that this 
channel for produ ci ng stellar sph e roids i s sign if icant (e.g.. 



Abadi et al.l I2OO6I: IZolotov et al.1 lioogl. l2010l: 



Oser et all 



2OIOI : IScannapieco et al.1 12OIII : iFont et al.1 l2011l . hereafter 
Paper I). An important caveat to bear in mind, is that 
many cosmological simulations suffer from a large degree 
of overcooling, and the importance of the in situ compo- 
nent could thus be overestimated in suc h simulations . How - 
ever, as we showed in Paper I and in ICrain et al.1 (j2010l ) 
hereafter CIO), th e GIMIC suite of cosmological simulations 
Grain et al] |2009| ) does not suffer from a large degree of 
overcooling on the scale of Milky Way mass (and lower) 
systems, due to efficient (but energetically-feasible) super- 
nova (SN) feedback. Note that strong SN feedback is widely 
believed to be required t o explain the generally low star 
formation efficiency (e.g ., iGroton et al.l I2OO6I : I Bower et al.l 
I2OO6I : Schave et al. 2010l ) and the 'missing baryon problem' 
(e.g.. 'Bre gmanll2007l : CIO) of normal galaxies, as well as the 
enrichm ent of the intergalactic medium (e.g., Aguirre et al.l 
[200 ll : iQppenheimer fc DavellioO^ : IWiersma et al.ll201ll ). As 



we demonstrated in Paper I, the gimic simulations yield 
satellite luminosity functions and stellar spheroid surface 
brightness distributions that are quite comparable to those 
observed for the Milky Way and M31, which is a direct result 
of efficient SN feedback. 

We also demonstrated in Paper I that the GiMiC simu- 
lations successfully reproduce the radial distribution of met- 
als (i.e., the large-scale gradient) in the spheroid, a problem 
that has been a thorn in the side of hybrid methods for some 
time. This begs the question of whether or not a significant 
in situ component can also alleviate the problems that hy- 
brid methods have in reproducing the shape of the spheroid, 
the evidence for a 'dual halo' from kinematics, and also the 
degree of substructure. 

In the present study, we focus on the global structure 
and kinematics of spheroids and leave the study of substruc- 
ture for a future work. In particular, in the present study we 
first examine in some detail the nature of the in situ compo- 
nent in the GiMiC simulations and show that its origins are 
deeply rooted in the dynamical heating of the proto-disc at 
high redshift. The implication of this finding is that the in 
situ component of the spheroid is expected to have differ- 
ent kinematic characteristics (namely rotation and perhaps 
also a different velocity anisotropy profile) than that of the 
spheroid component built from the disruption of satellites. 
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In addition, one expects a more flattened distribution for 
this in situ component, due to its net rotation. We make 
comparisons witfi measurements of the shape and kinemat- 
ics of the Milky Way, M31, and more distant galaxies. The 
upshot of this comparison is that the simulated disc galaxies 
have spheroids of the correct shape (oblate with a median 
axis ratio of « 0.6 at r < 30 kpc) and that the kinemat- 
ics show evidence for two components (due to in situ vs. 
accreted), as observed. 

The present paper is organised as follows. In Section [2] 
we describe the simulations and the selection of the sample 
of Milky Way-mass disc galaxies. In Section 3 we examine 
the origin of the in situ component of the spheroid in the 
simulations. In Section 4 we present our main results on 
the global structure and kinematics of the spheroids of the 
simulated galaxy population and make comparisons with ob- 
servations. Finally, in Section 5 we summarise and discuss 
our findings. 



2 SIMULATIONS AND SAMPLE SELECTION 

We use the Galaxies-Intergalactic Medium Interaction Cal- 
culation (gimic) suite of cosmological hydrodynamical simu- 
lations, which were carrie d out by the Virgo C onsortium and 
are described in detail bv lCrain et al.l (|2009l ). hereafter C09 
(see also CIO and Paper I). We present only a very brief 
summary of the simulations here and refer the interested 
reader to the above papers for further details. 

The suite consists of a set of hydrodynamical re- 
simulations of five nearly spherical regions (~ 20h~^ 
Mpc in radius) extr acted from the Millennium Simulation 
l|Springel et aLboOSl V The regions were picked to have over- 
densities at z = 1.5 of (-1-2, +1,0, —1, — 2)cr, where a is the 
root-mean-square deviation from the mean on this spatial 
scale. The 5 spheres are therefore environmentally diverse 
in terms of the large-scale structure that is present within 
them. For the purposes of the present study, however, we will 
only select systems with total 'main halo' (i.e., the dominant 
subhalo in a friends-of-friends, hereafter FoF, group) masses 
similar to that of the Milky Way, irrespective of which GiMiC 
volume it is in (i.e., irrespective of the large-scale environ- 
ment). C09 found that the properties of systems of fixed 
main halo mass do not depend significantly on the large- 
scale environment (see, e.g., Fig. 8 of that study). 

The cosmological parameters adopted for GiMiC are the 
same as those assumed in the Millennium Simulation and 
correspond to a ACDM model with Qm = 0.25, = 0.75, 
Qb = 0.045, as = 0.9 (where erg is the rms amplitude of 
linearly evolved mass fluctuations on a scale of 8h~^ Mpc 
at z = 0), Ho = WOh km s"^ Mpc"\ h = 0.73, = 1 
(where is the spectral index of the primordial power 
spectrum). The value of erg is approximately 2-sigma higher 
than has been inferr ed from the most recent CMB data 
(|Komatsu et al.ll2009l l. which will affect the abundance of 
Milky Way-mass systems somewhat, but should not signifl- 
cantly affect their individual properties. 

The simulations were evolved from z = 127 to z — 
using the Tree PM-SPH code GADGET-3 (last described in 
ISpringel|[2005l ). which has been modified to incorporate 
new baryonic physics. Radiative cooling rates for the gas 
are computed on an element-by-element basis by interpolat- 



ing within pre-comp uted tables generated with CLOUDY 
(|Ferland et al.lll998l ) that contain cooling rates as a func- 
tion of density, temperature, and redshift and that ac- 
count for the presence of the cos mic microwave background 
and photoionisation from a iHaardt fc Madaul (i20 01 .) ion - 
ising UV/X-Ray background fsee IWiersma et al.l l2009al V 
This background is switched on at z = 9 in the simula- 
tion, where it quickly 'reionises' the entire simulation vol- 
ume. Star formatio n is tracked in the si mulations following 
the prescription of iSchave fc Dalla Vecchia (20 0^) which re- 
prod uces the observed Kennicutt-Schmidt law (|Kennicuttl 
Il998l ) by construction. The timed release of individual el- 
ements by both massive (Type II SNe and stellar winds) 
and intermediate mass stars (Type la SNe and asymptotic 
giant branch stars) is i ncluded following the prescription of 
|w lersma et al.l (|2009bh . A set of 11 individual elements are 
followed in these simulations (H, He, C, Ca, N, O, Ne, Mg, 
S, Si, Fe), which represent all the important species for com- 
puting radiative cooling rates. 

Feedback from SNe is incorporated usin g the kinetic 
wind model of iDalla Vecchia fc Schav3 (|2008l ) with the ini- 
tial wind velocity, set to 600 km/s and the mass-loading 
parameter (i.e., the ratio of the mass of gas given a velocity 
kick to that turned into newly formed star particles), 77, set 
to 4. This corresponds to using app roximately 80% of the 
total energy avaUable from SNe for a lChabriej (^003) IMF, 
which is assumed in the simulation. This choice of parame- 
ters yields a good match to the peak of the star formation 
rate history of the universe (C09; see also lSchave et al1l2010l ) 
and reproduces a number of X-ray /optical scaling relations 
for normal disc galaxies (CIO). As shown in Paper I, the 
GIMIC simulations also produce realistic spheroidal compo- 
nents around ~ L, disc galaxies (see Figs. 3-5 of that paper, 
which compare the predicted surface brightness and metal- 
licity profiles to observations). 

To test for numerical convergence, the GiMiC simula- 
tions have been run at three levels of numerical resolution: 
'low', 'intermediate', and 'high'. The low-resolution simula- 
tions have the same mass resolution as the Millennium Sim- 
ulation (which, however, contains only dark matter) while 
the intermediate- and high-resolution simulations have 8 and 
64 times better mass resolution, respectively. Only the — 2cr 
volume was run to 2 = at high resolution, owing to the 
computational expense of running such large volumes at 
this resolution. As in Paper I, we adopt the intermediate- 
resolution simulations in the main analysis and reserve the 
high-resolution —2a simulation to test the numerical conver- 
gence of our results. The intermediate-resolution runs have 
a dark matter particle mass moM — 5.30 X lO'^h-^ Mq and 
an initial gas particle mass of mgaa ^ 1.16 X lO'^/i"^ M0, 
implying that it is possible to resolve satellites with total 
stellar masses similar to those of the classical dwarf galeixies 
around the Milky Way (with M* 10'^"^°Mq) with several 
hundred up to ^ 1000 particles. The gravitational soften- 
ing of the intermediate-resolution runs is lh~^ kpc. As we 
showed in the Appendix of Paper I, the surface brightness 
and metallicity profiles are converged at this resolution, as is 
the satellite luminosity function for satellites brighter than 
Mv ~ —12.5 (note that bright satellites dominate the as- 
sembly o f the accreted component of the stel lar spheroid; 
see, e.g.. lFont eraLll2006al : ICooper et al.ll2010l ). In the Ap- 
pendix of the present paper, we show the global shape and 
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kinematics of the spheroid in the intermediate-resolution 
runs agrees well with that of the high-resolution run. 

Note that for all these simulations the remainder of the 
(500/i~^ Mpc)"^ Millennium volume is also simulated, but 
with dark matter only and at much lower resolution. This en- 
sures that the influence of the surrounding large-scale struc- 
ture is accurately accounted for. 

2.1 Sample selection 

We use the same sample of 412 Milky Way-mass disc galaxies 
introduced in Paper I. A Milky Way-mass system is a de- 
fined as a system with a present-day total (gas+stars-|-dark 
matter) mass within r2oo (i-e., the radius which encloses a 
mean density of 200 times the critical density of the uni- 
verse) of 7 X 10" Af© < Ahoo < 3 X W^^Mq. This range 
roughly spans the mass estimates in the literature for both 
the Milky Way and M31, although we are not attempting to 
reproduce the properties of either of these systems in detail. 
Systems and their substructures are id entified in the simu - 
lations using the Subfind algorithm of [Oola^^^JiSSlL, 
that e xtends the standard implementation of Sprin gel et al.) 
l|200ll ) by including baryonic particles in the identification of 
self-bound substructures. For each FoF system that is iden- 
tified, a spherical overdensity (SO) mass with A — 200 is 
computed (i.e., M2oo)- The gas-f stars+dark matter associ- 
ated with the most massive subhalo of a FoF system are con- 
sidered to belong to the galaxy, while the gas-|-stars-|-dark 
matter associated with the other self-gravitating substruc- 
tures are classified as belonging to satellite galaxies. As the 
current study is focused on examining the extended stellar 
distributions of normal disc galaxies, we exclude all bound 
substructures (i.e., satellites) from the analyses presented in 
this paper. We also exclude mass that is part of the FoF sys- 
tem but which is not bound to it (the so-called FoF "fuzz"), 
but note that these particles make a negligible contribution 
to the stellar spheroid by mass. 

Galaxies identified in this way are then classified mor- 
phologically as disc- or spheroid-dominated, based on their 
dynamics. In particular, we use the ratio of disc stellar mass 
to total stellar mass (D/T) within 20 kpc computed by CIO 
to select our sample. In particular, we select disc galaxies 
with D/T > 0.3. We note that CIO decomposed the spheroid 
component from the disc by first computing the angular mo- 
mentum vector of all the stars within 20 kpc, calculating the 
mass of stars that are counter-rotating, and doubling it. This 
procedure therefore assumes that the spheroid has no net an- 
gular momentum. The remaining mass of stars is assigned 
to the disc and D/T is computed as the disc mass divided 
by the sum of the disc and spheroid masses. In reality, how- 
ever, the spheroid of the simulated galaxies (and likely of 
real galaxies as well) does in general possess a net angular 
momentum. We explicitly show this below for the simulated 
galaxies. However, ignoring this fact in the selection of our 
sample has no important consequences for our results or 
conclusions, since the vast majority of our simulated galax- 
ies are highly disc-dominated regardless of whether or not 
we take into account the angular momentum of the spheroid 
when calculating D/T. 

The method of CIO for decomposing the disc and 
spheroid does not assign individual star particles to either 
component, it simply computes the total mass in two com- 



ponents. We follow the approach of lAbadi et all (|2003h m 
assigning particles to the disc and spheroid. In particular, 
we align the total angular momentum of the stellar disc with 
the Z-axis, calculate the angular momentum of each star par- 
ticle in the X-Y plane, Jz, and compare it with the angular 
momentum, Jcirc, of a star particle on a co-rotating circular 
orbit with the same energy. Disc particles are then sel ected 
using a cut of Jz/ Jdic > 0.8 as in IZolotov et all (|2009l ). Vi- 
sually, this cut is quite successful in isolating the disc, but in 
a number of cases we found that particles well above/below 
the disc mid-plane were also assigned to the disc. For this 
reason, we also apply a spatial cut such that disc particles 
cannot exist more than 2h~^ kpc above or below the disc 
mid-plane (in Fig. [T] below we show that star- forming gas in 
the disc is generally confined to this region). We have also 
tried assigning disc particles based on the fracti on of stellar 
kinetic energy in ordered rotation, as defined in 'Sal es et al.l 
(2010) (sec also Fig. [2] below) , but found no significant dif- 
ferences from the method outlined above. 

Finally, we point out that we have not imposed any 
explicit constraints on the merger histories of the simu- 
lated galaxies. This is in contrast with many of the pre- 
vious cosmological modelling studies of Milky Way-type 
galaxies, which usually select their simulation candidates 
from systems that have not undergone any major mergers 
in the recent past (e.g., ICooper et al.ll2010l ). As noted in 
Paper I, the constraints adopted in previous hybrid stud- 
ies are overly conservative, given that the majority of the 
Milky Way-mass systems in GIMIC are disc-dominated sys- 
tems (see CIO). The resilience of galactic discs to mergers 
likely stems from the fact that the galaxies contain a sig- 
nificant dissipational gaseous component, w hich is able to 
maintain (or re-establish ) the disc (see, e.g.. iHopkins et al] 
l2009l :[ Moster et al.ll201ll ). In any case, for the present study, 
which is based on hydrodynamic simulations of large vol- 
umes, there is no need to place explicit constraints on the 
merger histories since it is easy to identify which galaxies 
have prominent discs at z — 0. 

Applying the system mass and morphology criteria de- 
scribed above yields a total sample of 412 systems in the 5 
intermediate-resolution GiMiC volumes. 

2.2 Distinguishing accretion and in situ formation 

We define in situ star formation as star formation that takes 
place in the most massive subhalo of the most massive pro- 
genitor (MMP) FoF group of the present-day system. Stars 
are said to be 'accreted' if they formed either in another 
FoF group or in a satellite galaxy (i.e., not the dominant 
subhalo) of the MMP. In practice, only a small proportion 
(2%) of the stellar halo is from stars that formed in satellites 
of the MMP and were then stripped; most of the accreted 
stars were formed in other FoF groups (i.e., other haloes) 
before they became satellites of the MMP. 

To identify the MMP of a given Milky Way-mass system 
at 2 = 0, we first select all dark matter particles within r2oo 
belonging to the dominant subhalo. We then use the unique 
IDs assigned to these particles to identify the FoF group in 
previous snapshots that contains the largest number of these 
particles. This earlier FoF group is defined to be the MMP 
and the dominant subhalo of that FoF group is assumed to 
be the most massive progenitor of the dominant subhalo of 
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Figure 1. Stacked spatial distributions (using all 412 galaxies) of the cold, star forming gas at various redshifts z^hs- Only gas in the 
most massive progenitor subhalo is shown (i.e., in situ star forming gas). In the top row the individual galaxies have been aligned before 
stacking so that the angular momentum vector of their stars at z = is parallel to the Z-axis. The median misalignment angle between 
the cold gas ai z = z^hs a^nd the stars at z = is quoted in the legend. In the middle and bottom rows the individual galaxies are aligned 
using the z = Zq^s angular momentum vector of the stars and cold gas, respectively, before stacking. The red and blue filled contour 
enclose 50% and 90% of the star forming gas particles, respectively. The vast majority of in situ stars formed in a disc. The median 
misalignment between the orientation of the disc at high redshift (where the majority of the in situ stars in the present-day spheroid 
formed, see Fig. [Sj and the disc today is large, implying that most discs have been severely torqued between z r^ 2 and the present day. 



the system at z = 0. Each of the star particles that comprise 
the 2 = stellar spheroid is tracked back in time to the FoF 
group and subhalo it belonged to, if any, when it formed 
from a gas particle. If, at the time of star formation, the 
star particle is a member of the dominant subhalo of the 
MMP, the star particle is said to have formed in situ. If it 
formed in another FoF group or in a non-dominant subhalo 
of the MMP, then it is said to have been accreted. 

As discussed in Paper I, there are sufficient snapshots 
output for each of the GIMIC volumes to be able to reliably 
assign particles to a formation mode. We trace each of the 
2 = stellar halo particles back in time as far as 2 = 10 to 
determine whether they formed in situ or were accreted. 



3 ORIGIN OF THE IN SITU COMPONENT 

Before examining the shape and kinematics of the stellar 
spheroid, we begin by examining in some detail the origin 
of the in situ component, which dominates the spheroid by 
mass in our simulations. We show below that there is a close 
physical connection between the in situ component of the 
spheroid and the early disc. This motivates us to look at 
potential differences in the shape and kinematics of the in 
situ and accreted components of the spheroid in Section 4. 



In Fig. [T] we show the stacked (using all 412 disc galax- 
ies) spatial distribution of cold, dense star forming gas 
(T < 10^ K and nu > 0.1 cm""') at various redshifts, 2obs. 
We only include gas that is part of the most massive sub- 
halo of the most massive progenitor of the disc galaxies that 
were selected at z = 0. This is, by definition, the gas that 
gives rise to the in situ stellar component of the galaxy (both 
of the disc and the spheroid). The red and blue filled con- 
tour enclose 50% and 90% of the star forming gas particles, 
respectively. 

In the top row we show the stacked distributions where 
the individual galaxies have been aligned before stacking so 
that the angular momentum vector of their stars at 2 = 
is parallel to the Z-axis. At 2 = (top left panel) a disky 
distribution is evident. This implies that the bulk of the star 
formation is occurring in a disc and that the star forming 
gas disc is well-aligned with the stellar disc. As we track 
the most massive progenitor of these galaxies back to pre- 
vious redshifts and stack using the same (z = 0) rotation 
matrix to align the gas the stacked star forming region be- 
comes less disky and more isotropic in appearance. Either 
the star forming gas was genuinely distributed in this way 
(and therefore discs only formed very recently), or the gas 
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was in a disc but the disc has been significantly torqued 
between z — and z — Zobs- 

We can distinguish between these two possibilities by 
aligning the star forming gas at Zobs using either the an- 
gular momentum of the stars (middle row) or the angular 
momentum of the cold gas itself (bottom row) a.t z — Zobs- 
A disky distribution for the star forming gas is readily ap- 
parent using either approach. The radial extent of the disc 
becomes smaller at high redshift which is qualitatively con- 
sistent with previous simulati on studies of disc galaxies as 
well as with observations fe.g.. lWang et al1l201ll ). Thus, the 
star forming gas becomes more isotropic looking with in- 
creasing redshift in the top row because of the significant 
torquing of the disc. In the legend in the top row we quote 
the median misalignment angle (in deg.), ^misalign, between 
the angular momentum vector of the star forming gas at 
z = Zobs and that of the stars at z — 0. A random distribu- 
tion (i.e., no correlation between the vectors) would have a 
median misalignment of 90 degrees. This is nearly the case 
of the star forming discs at 2 = 2.0, implying that there is al- 
most no correlation between the orientation of (proto)-discs 
at 2 = 2.0 and the descendant discs today. 

Based on Fig. [1] we therefore conclude that the vast 
majority of the overall (disc+spheroid) in situ component 
present in our galaxies at 2 = formed in a disc at some 
previous time. Since stars formed in situ dominate the stellar 
spheroid in these simulations out to rso ~ 20 — 30 kpc (see 
Fig. 7 of Paper I), this implies that some process has added 
sufficient energy ('puffed up') to the stars to liberate them 
from the confines of the disc. We can explicitly check that 
this is the case by comparing the ratio of (specific) kinetic 
energy in ordered rotation [Krot = {Jz / R)'^ /2, where R = 
(X^ -I- y2-)i/2j ^^^^Y kinetic energy [Ktot = where 
vtot = {vx +Vy + '^1)^''^] for in situ stars at the present-day 
to the same ratio for the same particles at the time of star 
formation. Operationally, we select cold, star forming (in 
situ) gas particles residing in the most massive progenitor at 
some previous redshift Zobs and then find the corresponding 
star particles at 2 = using the unique particle IDs stored in 
the snapshot data. Note that not all of the cold, star forming 
gas turns into stars by 2 = (e.g., some gas is heated up 
by SNe) and not all of the in situ stars at 2 = will be in 
the cold gas component at 2obs. We are therefore examining 
a subset of the particles, but it is an unbiased selection in 
terms of kinematics. 

The results of this test are shown in Fig. [S] At any 
redshift 2obs the cold, star forming gas has a much larger 
fraction of its energy in highly ordered rotation (e.g., 
Krot/Ktot > 0.8) than the stars at 2 = which formed 
from that gas. As we look at cold gas further back in time 
(increasing 2obs), we see that the stars at 2 = that formed 
from that gas have lower and lower fractions of their energy 
in highly ordered rotation. This is what one would expect, 
as there has been more time for the stars to be dynamically 
puffed up (e.g., via accretion, minor mergers, disc instabili- 
ties, and so on). 

When did the stars in the in situ spheroid form? In 
Fig. Owe show the stacked distribution of stellar ages of the 
disc and the in situ and accreted spheroid components at 
2 = 0. As expected, the accreted spheroid is old (median 
age of Si 11.1 Gyr, or z ^ 2.7), having been built up from 
the tidal disruption of dwarf galaxies. As shown in Paper 



stars a: z = 

COS at z = Zq(,5 ' 

= 0.5 ; 


: =1.0 1" 

: 

1 r 


1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 

Zobs =1-5 1 


1 1 1 1 1 1 1 1 1 r 1 1 1 1 1 1 1 1 1 

; z„,3 = 2.0 1" 



0.0 0.? 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0 

K.ot/K.ot 



Figure 2. The fraction of kinetic energy in ordered rotation for 
in situ stars at z = and the cold gas which formed them at 
z = Zobs- The different panels show different values of Zoba- The 
specific kinetic energy in ordered rotation is Krot = {Jz / R)'^ 
(where R = (X'^ + Y'^)^^'^) and the total specific kinetic energy is 
Ktot = -"tot/^ (where i;tot = (vj^ + v'^ + v^^)^/'^). Operationally, 
cold, star forming (in situ) gas particles are selected at Zobs and 
the corresponding star particles at z = are found using the 
unique particle IDs. At any redshift Zobs the cold, star forming 
gas has a much larger fraction of its energy in highly ordered 
rotation (e.g., i^rot/^S'tot > 0.8) than the stars at z = which 
formed from that gas. This difference increases with increasing 
redshift. 

I, however, the accreted spheroid is still being assembled 
today (see Fig. 10 of that study). Both the in situ spheroid 
and disc are on average younger than the accreted spheroid. 
Interestingly, the median age of the in situ spheroid (« 8.0 
Gyr, corresponding to z « 1.1) is significantly older than 
the median age of stars in the present-day disc (« 6.2 Gyr, 
or 2 « 0.7). Taken together with Figs. 1-2, this implies that 
the in situ component of the spheroid primarily formed in a 
proto-disc that evidently was heated up before much of the 
current disc stellar population was formed. 

Interestingly, there is a marked decline in the distri- 
bution of the in situ spheroid toward younger ages, even 
though the distribution of ages for the disc is effectively flat 
below « 9-10 Gyr (implying a nearly constant disc star for- 
mation rate for the past « 9-10 Gyr). This suggests that 
the heating mechanism that is responsible for producing the 
in situ spheroid is becoming less and less efficient at late 
times. A likely culprit for the heating is the accretion of 
satellites (an d mass in general) onto the simulated galax- 
ies. Recently, iBovlan-Kolchin et al.l l|2010l ) studied the mass 
accretion histories of thousands of Milky Way-mass dark 
matter haloes in the Millennium-II simulation. If one de- 
fines the epoch of formation as the time whe n the halo has 
accre t ed half of its present-day m ass (e.g., iLacev fc Cold 
Il993l ). [Bovlan-Kolchin et al.l l|2010l ) find a median epoch of 
formation of z « 1.2 for Milky Way-mass haloes. This is 
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Figure 3. Distribution of stellar ages. Most of the accreted com- 
ponent of the spheroid is very old (typically 11-12 Gyr). The 
disc is younger and has roughly had a constant SFR for the past 
9 — 10 Gyr. The in situ spheroid component has a slightly older 
median age than the disc. The SFR of this component has been 
declining steadily since z 1.5. Note that the amplitudes of the 
histograms are determined by normalising to the total number of 
stars in each population separately. 

suggestively close to the median age of the in situ spheroid 
component. Alternatively, instead of using the epoch of for- 
mation, one can use the re dshift of most massive merger. 
iBovlan-Kolchin et all (|2010l ) find a median mass ratio of 10:1 
for the most massive merger (i.e., the mass of the satellite 
at accretion is 10% of the final virial mass of the main halo). 
We have tracked in the GIMIC simulations all of the galaxies 
in our sample which have undergone a merger with mass 
ratio 10:1 (or lower). The median redshift for systems with 
such mergers is ~ 1.1. Both arguments point to dynamically 
active systems (on average) at this period of tim^fl- After 
z ^ 1 th e mass accretion histories fi atten significantly (see 
Fig. 2 of lBovlan-Kolchin et aLlbOloD . which is qualitatively 
consistent with the decline in the stellar age distribution of 
the in situ spheroid. 

Our analysis th erefore favours a scen ario similar to that 
proposed recently bv lPurcell et~aLl (|2010l ). in which the inner 
parts of galactic spheroids contain proto-disc stars that may 
have been liberated in the same merger/accretion events 
that delivere d mat erial to the outer spheroid. Note that 
iPurcell et al.l l|2010l ) primarily analysed idealised collision- 
less simulations of disc bombardment by infalling satellites, 
whereas GiMiC is a suite of fully cosmological hydrodynamic 
simulations. 

Finally, we make a brief comment on the thermody- 

* The lack of a correlation between the orientation of the present- 
day disc and proto-disc in which the in situ spheroid stars were 
formed (evident in Fig. [TJ adds further evidence of a dynamically 
active state at z ~ 1. 



Figure 4. Distributions of maximum past temperatures of the 
gas which formed the disc and the accreted and in situ spheroid 
components. The black, blue, and red dotted vertical lines corre- 
spond to the median Tmax for the disc and the accreted and in situ 
spheroid components, respectively. The accreted spheroid formed 
mostly from cold gas ( "cold mode" ) , whereas approximately half 
of the stars in the disc and half in the in situ spheroid formed 
from gas that was accreted in the so-called "hot mode" (i.e., they 
formed from the cooling of gas that was first shock heated to 
approximately the virial temperature) . 

namic state of the gas from which the in situ component 
formed. In Fig. [3] we show the distribution of the maxi- 
mum past temperature of the gas particles from which the 
disc, in situ spheroid and accreted spheroid star particles 
formed. The distributions are obviously bimodal. The dis- 
tributions for the disc and in situ spheroid components are 
quite similar, with a pealjf] at Tmax ~ 10^ K and another 
at Tmax ~ ]^q6.o-6.3 Corresponding to so-calle d "cold 
mode" and "hot mode " accretion (e.g.. [Blrn boim fc Dekell 
l2003l : iKeres et al.ll2005l . l2009l : Ivan de VoorTet al.ll2oTlh . re- 
spectively (note that the virial temperature for haloes we are 
studying here is Tvir ~ 10^ K). In agreement with CIO, we 
find that approximately half of the disc stellar mass formed 
via the cooling of hot gas. We find a similar result for the 
in situ spheroid component. By contrast, the accreted com- 
ponent formed predominantly from gas that has never been 
hot. 

The above points to a different physical n ature for the in 
situ co mponent than advocated previously bv lZolotov et alj 
(|2009l ). In particular, these authors found a median red- 

^ This lower Tmax peak corresponds to the temperature floor im- 
posed by the photoionising UV/X-ray background in the simula- 
tions. Note that the temperature floor that is established depends 
on redshift, as the cooling rate of the gas increases with redshift 
due to the increase in density with redshift. This is why the lower 
Tmax peak for the accreted co mponent, which formed a t higher 
redshift, is only ^ lO*-^ K (see lvan de Voort et al]|201ll) . 
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shift of formation of ~ 3 for the in situ spheroid compo- 
nent in their simulations. Furthermore, the stars formed at 
very small radii (< 1 kpc) from gas that was fed via cold 
mode accretion. In contrast, we find the in situ component 
is younger, formed in a proto-disc, and that approximately 
half of the stellar mass was formed via the cooling of gas that 
was heated to approximately the virial temperature of the 
system. A possible reas on behind this differe nce is that the 
simulations adopted by IZolotov et al .' ( 2009') suffered from 
a larger degree of overcooling due to inefficient feedback. A 
likely consequence of this overcooling is that a significant 
fraction of the gas that would have been destined to end 
up as part of the in situ stellar halo was instead converted 
into stars in small dwarf galaxies at high redshift, and were 
therefore instead incorporated into the accreted spheroid 
component when these systems fell into the main progenitor. 
This sce nario is consist ent with the fact that the simulations 
used by IZolotov et al.l (,2009) produce present-day satellite 
galaxies that are too brig ht compared to those of the Milky 
Way or M31 (see Fig. 4 of lGovernato et al.llioOTi '). We note, 
however, that there are many other differences between the 
simulations that could also play a role (e.g., somewhat dif- 
ferent implementations of star formation, chemodynamics, 
and metal-dependent radiative cooling). Direct comparisons 
of the properties of stellar spheroids predicted by difi'erent 
simulation codes all using the same initial conditions would 
obviously be beneficial in helping to firmly identify the ori- 
gin of the apparent differences between the simulations. 

The above results point to a close physical connection 
between the disc and the in situ spheroid in our simulations. 
Naively, we therefore expect that the latter should have kine- 
matics and a spatial configuration that may be more remi- 
niscent of the former than to the accreted spheroid. As we 
highlighted in Section 1, there is mounting observational ev- 
idence for net rotation and flattening of the inner spheroid of 
both the Milky Way and M31. This is a tantalising sugges- 
tion that the inner halo indeed contains a significant dissipa- 
tional component. Below we study the shape and kinematics 
of the stellar 'spheroid' and, where possible, make compar- 
isons to the observational data. 



4 STRUCTURE AND ROTATION 
4.1 Structure 

Below we examine the structure of the spheroid from exter- 
nal and internal perspectives separately. We make compar- 
isons to measurements of M31 and NGC 891 and stacked 
observations of more distant galaxies when considering the 
former and to the Milky Way when considering the latter. 



4-1.1 Structure - External view 

We begin by examining the average 2D stellar distribution of 
the spheroid. To this end, we stack all 412 simulated galax- 
ies in an edge-on configuration, by first aligning the angular 
momentum vector of the stars of each galaxy with the Z- 
axis of the simulation box. Since we are interested in the 
diffuse spheroid, we stack only the star particles that are 
bound to the most massive subhalo (the main galaxy; i.e.. 



we exclude star particles bound to orbiting satellite galax- 
ies). To facilitate easy comparisons with observations (see 
below), we compute a stacked SDSS i-band surface bright- 
ness matF°l. The particles are interpolated to a 2D regular 
grid after stacking using the smooth particle hydrodynamics 
(SPH) kernel with a large number of smoothing neighbours 
(512). This was done to smooth out substructures not asso- 
ciated with self-bound satellites (again we are interested in 
the global/average properties). We have checked that using 
smaller numbers of smoothing neighbours does not aff'ect 
the average shape as a function of projected radius, it does 
however result in a more noisy profile. 

The resulting smoothed i-band map is shown in the left 
panel of Fig. [5] The numbers associated with the contours 
give the mean surface brightness in mag arcsec"^. A zoom 
in on the central 50 kpc is shown in the right panel. Note 
that the disc stars are located at \Z\ < 2 kpc (see Fig. [l|. 
It is readily apparent that the spheroid component shows 
a significant degree of fiattening and that the fiattening is 
muc h more pronounced in the inner regions. 

IZibetti et all \2004 ) have performed a stacking analy- 
sis of approximately 1000 edge-on disc galaxies in the SDSS 
and have measured the surface brightness distribution of 
the spheroid down to ~ 30 mag arcsec"^ in the i-band. 
We can directly co mpare our results above to this study. 
IZibetti et al.l (|2004l ) have quantified the degree of fiatten- 
ing by fitting ellipses to their stacked data at a variety of 
surface brightness levels. Their trend oi b/a (where b and a 
are the semi-minor and semi-major axes respectively) with 
distance along the major axi^^ is shown as the solid red 
curve in Fig. [G] We have performed an identical analysis 
to the stacked i-band image of all 412 simulated galaxies, 
which is represented by the solid black curve. The dashed 
black curve represents the average fiattening profile when 
we stack only gal axies that would make the sample selec- 
tion criteria of Zi betti et al. ||2004), namely that the isopho- 
tal semi-major axis a exceed 10 arcseconds and that the 
isophotal b/a be sm all in order to select disc galaxiej^. 
IZibetti et all HqoI) adopt an isophotal b/a threshold of 
< 0.25, whereas we adopt a slightly larger threshold of 0.4 
in order to boost our sample of galaxies. Only 41 galaxies 
in our sample have b/a < 0.25, whereas approximately half 
(210) have b/a < 0.4. Imposing a cut of b/a < 0.25 gives re- 
sults that are consistent with (but considerably more noisy 
than) our default cut of b/a < 0.4. 

Overall there is good agreement in the shape of the 
fiattening profile between the simulations and the measure- 



" The i-band luminosities are computed by treating each star 
particle as a simple stellar population (SSP). The simulations 
adopt a Chabrier IMF and store the age and metallicity of the 
particles. We use this information to compute a spectral energy 
di stribution for each s tar p article using the GALAXEV model 
of iBruzual fc Charloti ||2003D . The i-band luminosity is obtained 
by integrating the product of the SED with the SDSS i-band 
transmission filter function. Our calculation neglects the effects 
of dust attenuation, which may be relevant at very small radii. 

We have converted their distances from pixel lengths into kpc 
by adopting the mean redshift of their sample, which is z = 0.05. 
^■^ Isophotal values of a and b refer to the values produced by the 
SDSS reduction pipeline at a surface brightness level of 25.0 mag 
arcsec"^. 
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Figure 5. Stacked SDSS i-band surface brightness contour maps of tlie simulated galaxies. Tiie galaxies are rotated to an edge-on 
configuration before stacking. Self-gravitating substructures (satellites) have been excised and the map has been heavily smoothed using 
an SPH kernel (see text for details). The numbers associated with the contours give the mean surface brightness in mag arcsec"^. Left: 
Large-scale view. Right: Zoom in on the central regions. Note the stellar disc is generally confined to |Z| < 2 kpc. It is readily apparent 
that the spheroid shows a significant degree of flattening, particularly in the inner regions (r < 40 kpc). 
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Figure 6. The flattening ib/a) of the average surface brightness 
distribution as a function of distance along the major axis (R). 
The solid red curve represents the flattening profllc derived by 
IZibetti et al.l 12004) from a stacking analysis of Ri 1000 edge-on 
disc galaxies in the SDSS. The solid black curve is the flattening 
profile derived from an identical stacking analysis of all 412 simu- 
lated disc galaxies. The dashed black curve represents the average 
flattening profile when we stack only galaxies that would make 
the sample selection criteria of IZibetti et al. 1 ll2004h (see text). 
The flattening proflle of the simulated galaxies bears a remark- 
able resemblance to that derived from observations. 



ments of lZibetti et al.l (|2004l ). This statement is independent 
of whether or not we a pply the same select ion criteria for the 
disc galaxy sample as lZibetti et al.l (|2004l ). When we apply 
the same selection criteria, however, this results in an im- 
proved match to the normalisation. At r < 10 kpc the profile 
asymptotes to a level b/a « 0.2. This is due to the presence 
of the stellar disc. Beyond this radius we begin to probe into 
the spheroid. Even out to radii of 50 kpc the spheroid is still 
significantly flattened (b/a ~ 0.7). 

The agreement between the simulations and observa- 
tions is non-trivial, as it has proved notoriously difficult to 
produce reasonable disc galaxies at all in cosmological hy- 
drodynamic simulations. As discussed in Paper I, part of 
the reason for the improvement with GiMiC is that it does 
not suffer from strong overcooling for Milky Way-mass (and 
lower) galaxies owing to efficient (but energetically feasible) 
SN feedback. The feedback parameters themselves have been 
tuned only to match the peak of the global star formation 
rate density of the universe. As noted in the Introduction, 
strong SN feedback is also required to explain the low baryon 
fractions of normal galaxies and the enrichment of the in- 
tergalactic medium. 

In the above we have focused on the average shape of the 
spheroid that would be derived by stacking large samples of 
galaxies. This was required to probe the shape of the outer 
spheroid (beyond r > 30 — 40 kpc), since our simulations 
contain too few particles to probe the shape to very large 
radii on a system-by-system basis. As we will show later, 
however, the simulations predict that there should be a sig- 
nificant degree of system-to-system scatter in both the shape 
and kinematics of the inner spheroid, where observations are 
generally currently limited to anyway. At present it is quite 
challenging to probe the spheroids of individual galaxies. 
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Figure 7. Stacked stellar mass density contours in cylindrical 
coordinates. Self-gravitating substructures (satellites) have been 
excised and the map has been heavily smoothed using an SPH 
kernel (see text for details). The contours reflect the logarithm 
of the stellar mass density in Mq kpc"'' (for the spherically- 
averaged stellar mass density refer, to Fig. 3 of Paper I). The 
various panels show the density distribution for different cuts in 
the metallicity and age of the star particles. A comparison of the 
top and bottom panels shows that the flattening of the stellar 
mass density is much more pronounced for younger and more 
metal-rich stars and at smaller radii, in qualitative agreemen t 
with the measurements of the Milky Way bv lCaroUo et al] ||2007|) . 



owing to the very low surface brightness of this component. 
While in the future this will be readily achievable for large 
numbers of galaxies with, e.g., the LSST, studies of indi- 
vidual galaxies a re presently limited to the l ocal v olume. In 
the case of M31, IPritchet fc van den Berghl (|l994l 'l find that 
the spheroid is significantly flattened with h/a ~ 0.55 at a 
minor axis distance of 10 kpc. More recently, llbata et al.l 
have found evidence for an extended disc-like struc- 
ture out to 50 kpc which is aligned with the stellar disc 
and has b/a ~ 0.6. The large panoramic view of M31 to 
be provided by the ongoing Pan- Andromeda Archaeological 
Survey (PAndAS) will be useful for quantifying the degree of 
flattening as a function of distance out to very large radii. 
The flattening of the spheroid of NGC 891, a loca l Milk y 
Way analogue, has been measured by llbata et al.l l|2009l ). 
They find a fiattening of b/a ~ 0.5 out to a distance of ~ 25 
kpc. Both M31 and NGC 891 therefore show a level of fiat- 
tening that is comparable to the mean value derived from 
our simulated galaxies. 



set to take another leap in the coming years with the re- 
lease of data from SDSS-III/SEGUE-2 (Rockosi et al. in 
prep), RAVE, SkyMapper, Pan-Starrs and the launch of 
Gaia. Some of the best constraints on the formation and 
evolution of spheroids and galaxies in general will therefore 
come (and are now coming) from our own Galaxy. With 
this as motivation, we have examined the structure of the 
spheroid from an internal perspective. 

In Fig. [7] we plot stacked stellar mass density con- 
tours in cylindrical coordinates. The maps again exclude 
self-gravitating substructures (satellites) and have been 
smoothed to a 2D regular grid (in R,Z space) after stacking 
using the SPH smoothing kernel with 512 smoothing neigh- 
bours. The contours reflect the logarithm of the stellar mass 
density in Mq kpc~^ (for the spherically-averaged stellar 
mass density, refer to Fig. 3 of Paper I). The various panels 
show the density distribution for different cuts in the metal- 
licity and age of the star particles [i.e., analog ous to maps 
made of the Milky Way bv lCarollo erall (|2007l ) using SDSS 
data, see their Fig. 5]. 

A comparison of the top and bottom panels reveals 
that the flattening of the stellar mass density is much more 
pronounced for younger and more metal-rich stars and at 
small er radii. Qual i tative ly this is in good agreement with 
what ICaroUo et all (|2007l ) flnd for the Milky Wa£f|. By se- 
lecting young, metal-rich stars, we are predominantly select- 
ing the stars which formed in situ in the simulated galax- 
ies (old, metal-poor stars, by contrast, were predominantly 
formed in dwarf galaxies which have subsequently been ac- 
creted by the main galaxy and tidally destroyed). The flat- 
tening of the overall distribution, as clearly seen in Fig. [S] is 
therefore driven primarily by the in situ component of the 
spheroid. The gradual transition from a highly flattened to 
more spherical conflguration with distance from the galaxy 
centre (as in Fig. (Sjl mainly reflects the increasing impor- 
tance of the accreted component of the spheroid with in- 
creasing radius (see Fig. 7 of Paper I). 

In the above we have considered only the mean structure 
of the spheroid, derived by stacking all of our simulated 
galaxies. How much system-to-system scatter is there in this 
structure? We can assess this by fltting a smooth parametric 
model to the stellar mass distribution of each of the galaxies 
and then examining the scatter in the best-flt parameters 
of that model. Here we adopt an oblate single powerlaw 
distribution of the form 



p{R,Z) = 



Po 



[i?2 + (Z/g)2]"/2 



(1) 



where q = c/a and c and a are parallel and perpendicular 
to the angular momentum vector (Z-axis), respectively. This 
distribution has been commonly ad opted for modell ing the 
stellar halo of the Milky Way (e.g.. |jurfc et al.l l200^. Note 
that it has recently been shown that this model does not pro- 
vide an a cceptable fit to the Milky Way's outer halo (r > 25 
kpc, e.g.. lSesar et al.ll201ll : lD"eason et al.ll2011bl ). in that the 



4- 1-2 Structure - Internal view 

Measurements of the spatial structure (and substructure) 
and kinematics of the spheroid of our own Milky Way galaxy 
have rapidly increased in quantity and precision in recent 
years, particularly with the advent of the SDSS, and are 



Note we did not att empt to adopt the exact same cuts in 
metallicity employed by ICaroUo et al.l ||2007| ) since the simula- 
tions adopt empirical yields and Type la SNe rates, both of whic h 
arc uncertain at the factor of 2 level jWiersmaet al.ll2009bh . 
Thus, the metallicities of star particles in the simulations can 
be shifted up or down by at least 0.3 dex. 
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Figure 8. Scatter plot of the best-fit parameters (flattening c/a 
and power-law index a) of the oblate power-law (eqn. 1) derived 
from fitting to each of the simulated galaxies. For comparison we 
show recent observational estimates for the Milky Way. There is 
a large spread in the parameters but with a clear correlation be- 
tween the two. The estimates for the Milky Way are fully compat- 
ible with those of the simulated galaxy population, unlike hybrid 
models which lack any dissipational component and find highly 
prolate (as opposed to oblate) distributions for the spheroid. 



observed halo appears to drop off more quickly beyond this 
radius. Our simulations lack the mass resolution to probe 
the detailed shape of the outer halo on a system-by-system 
basis. However, as we showed in Paper I, the stacked mass 
distribution of our simulated galaxies does indeed drop off 
faster beyond r « 25 kpc out to r « 40 kpc, in rough accor- 
dance with observations. This 'kink' in the simulations rep- 
resents the transition from the in situ-dominated spheroid 
to the accretion-dominated spheroid. Beyond r 40 kpc 
the mass density profile asymptotes to a logarithmic slope 
of PS —3.5. In any case, for the purposes of the present study, 
we fit oblate single power-laws to the mass distribution to as- 
sess the system-to-system scatter in the simulations. When 
comparing to observations, we compare to the best-fit pa- 
rameters derived by fitting the same parametric model to 
the observational data over the same (or similar) range of 
radii. 

When fitting equation (1) to the simulated galaxies, we 
exclude star particles within \Z\ < 4 kpc (to exclude disc 
particles) and beyond a 3D radius rsD > 40 kpc. This selec- 
tioi£3 is meant to crudely mimic the typical selection criteria 
for stars in the halo of the Milky Way. Finally, before fit- 



We find that the best-fit parameters are generally insensitive 
to the outer radius cut but are somewhat sensitive to the cut 
in Z if no other cut is applied to exclude disc stars. In analysis 
of observational data, cuts on metallicity and/or kinematics are 
usually adopted to exclude disc stars. Furthermore, observations, 
e.g. with the SDSS, are limited in any case to > 30 deg, where 



ting the parametric model, each galaxy is rotated so that 
the direction perpendicular to the disc plane lies along the 
Z-axis. Thus, c/a < 1 implies flattening in the same sense as 
the disc, whereas c/a > 1 (which is physically and mathe- 
matically possible) would imply the halo is elongated in the 
direction perpendicular to the disc plane. 

In Fig. |5]we show a scatter plot of the best-fit parame- 
ters c/a and a where each black circle represents a galaxy. 
(In Section 4.3 we will explore how these parameters corre- 
late with the rotation of the spheroid and its in situ mass 
fraction; see Fig. 1121 1 It is readily apparent that there is 
a large range in the values of both the fiattening c/a and 
the power-law index a and that the two are strongly corre- 
lated, in that the density drops off more slowly with radius 
in highly flattened systems. The fiattening c/a ranges from 
roughly 0.25 to 1.4 with a median of 0.61. The power-law 
index ranges from roughly 2 to 4 with a median of 2.85. (For 
reference, in Paper I we found that the spherically-averaged 
logarithmic slope ranged from roughly 2.5 to 3.4 over a sim- 
ilar range of radii.) For comparison with the Milky Way, we 
show the best-fit parameters derived by fitting e quation (1) 
to (i) main-sequence turn-off stars in the SDSS (|Jurfc et al.l 
I2OO8I): (ii) blue horiz ontal-branch (BHB) stars in the SDSS 
( Deason et al.ll2011bl): and (iii) ma in-sequence turn-off stars 
in the CFHTLS (lSeiareE^[20lJ). In spite of the different 
samples and analysis techniques, similar results are obtained 
for the spheroid of the Milky Way in all three studies and 
all are quite compatible with the distribution of spheroid 
shapes in the simulations. 



4-1.3 Summary of 

Our analysis of the shape of the spheroidal component of 
the simulated galaxies in GIMIC indicates that spheroids are 
oblate in shape (in the same sense as the disc) and that 
the axial ratios derived from the simulations are consistent 
with measurements of individual nearby galaxies as well as 
with stacked observations of more distant (edge-on) galaxies. 
This is in stark contrast to the results of hybrid models for 
the formation of stellar haloes, which lack a dissipational 
component (i.e., the haloes are built from tidal disruption of 
satellites alone) and typically find a high ly prolate shape for 
the stellar halo (e.g.. lCooper et al] ^2010). The success of the 
GIMIC simulations in reproducing the shapes is a direct result 
of both the large fraction of the spheroid that is produced in 
situ and quite likely the timing of its formation as well. We 
speculate that if the in situ component had formed at very 
high redshift, well before the formation of a proto-disc, the 
shape of the spheroid would not be oblate with the observed 
axis ratio. 



4.2 Rotation 

4-2.1 Rotation - External view 

In Fig. [9] we show a stacked line of sight velocity (i.e., 
stacked along the y-axis) contour map, where all 412 simu- 
lated galaxies have been rotated to an edge-on configuration 



is the angle between the plane of disc and the height above the 
plane. 
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Figure 9. A stacked line of siglit {vy) velocity contour map. All 412 simulated galaxies have been rotated to an edge-on configuration 
before stacking. Red contours indicate motion away from the observer, whereas blue contours indicate motion towards the observer. The 
numbers associated with the contours give the line of sight velocity in km/s. Left: Large-scale view. Right: Zoom in on the central regions. 
Note that the stellar disc is generally confined to |Z| < 2 kpc. The 'spheroid' shows a significant degree of rotatio n, particularly in the 
inner regions (r < 40 kpc) and bears a close qualitative resemblance to the extended stellar disc around M31 (e.g.. llbata et a^l2005^ . 



before stacking. As in previous plots, the star particles are 
interpolated to a 2D regular grid after stacking using the 
SPH kernel with 512 smoothing neighbours. Red contours 
indicate motion away from the observer, whereas blue con- 
tours indicate motion towards the observer. The numbers 
associated with the contours give the line of sight velocity 
in km/s. The left panel gives a large-scale view and the right 
panel zooms in on the central 50 kpc. 

From the colour coding of the contours alone it is evi- 
dent that, on average, the spheroid has a net prograde ro- 
tation (i.e., in the same sense as the disc). At very large 
radii (r > 30 kpc) the mean rotation velocity is only 10-20 
km/s. In the inner regions, however, the spheroid rotates 
much faster and has the appearance of an extended, 'puffed 
up' disc that bears a close resemblan ce to the 'vast ext ended 
stellar disc' around M31 reported by (jlbata et al.ll2005l ). The 
lag between the disc and the 'puffed up' disc depends on the 
distance along the major axis but typical values range from 
-50 to -100 km/s, which is similar to what is found for M31. 

In Fig.[TD]we show the distribution of the ratio of the ro- 
tation velocity to the 3D velocity dispersion of the spheroid 
of the simulated galaxies (i.e., we compute one value of 
Kot/o" for each simulated galaxy's spheroid) for different 
cuts applied to the spheroid. For comparison, in each panel 
the black histogram represents the system-to-system scatter 
in the Kot/o" for the stellar disc. Note that in computing a 
single value of Vrot /o" for each galaxy we are averaging over 
all of the star particles in that galaxy (subject to the cut 
applied) and each particle is weighted equally in the averag- 
ing. We are thus examining the system-to-system scatter in 
the global kinematics of stellar spheroids in Fig. 1101 

In the top left panel of Fig. [10] we compare accreted and 
in situ spheroid components. Typically, the accreted compo- 



nent shows only a mild degree of prograde rotation (though 
in some cases the accreted spheroid is counter-rotating with 
respect to the disc) and is usually dispersion dominated 
(Kot/o" << 1). The in situ spheroid, by contrast, always 
shows co-rotation with the disc and shows a large spread in 
Kot/o", indicating that in some cases this component is dis- 
persion supported and in other cases it is rotation supported 
(on average Kot/o" 1). 

The origin of the rotation of the in situ component 
is clear. But it is interesting that the accreted component 
of the spheroid also typically shows some degree of pro- 
grade rotation with respect to the disc. Why might this 
be the case? Cosmological simulations have shown that 
there is preference for accretion of satellites along the ma- 
jor axis of the dark matt e r halo (e.g., throug h filaments; 
see iLibeskind et al. I l2005l : IZentner et al] 12005! ). Recently, 
iDeason et al.l ( 201 id ) used the GIMIC hydrodynamic simula- 
tions to show that there is also a general tendency (although 
with large scatter) for the plane of the galactic disc to lie 
along the major axis of the dark matter halo. Therefore, 
there is some preference for satellites to be accreted parallel 
to the plane of the disc, which could induce a net rotation of 
the accreted spheroid in the same sense as the disc. Another 
contributing factor may be the fact that galaxies on pro- 
grade orbits are preferentially susceptible t o orbital decay 
due t o dynamical friction (as argued by, e.g. jMorrison et al.l 
iooi), although this is only expected to be important for 
those satellites which are able to penetrate very close to the 
disc before being disrupted. But note it has even been ar- 
gued that th is effect can give rise to a small net retrograde 
motion (e.g.. iQuinn fc Goodmanlll986l ). 

In the other panels we show the distribution of Kot/o" 
for the spheroid component when we apply cuts in metal- 
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Figure 10. The distribution of the ratio of the rotation velocity to the 3D velocity dispersion of the spheroid of the simulated galaxies 
(i.e., one value of Vrot/ry is computed for each simulated galaxy) for different cuts applied to the spheroid. Top left: For the spheroid 
decomposed into in situ and accreted subcomponents. Top right: For the spheroid decomposed into relatively metal-rich and metal-poor 
subcomponents. Bottom left: For the spheroid decomposed into relatively young and old subcomponents. Bottom right: For the spheroid 
decomposed into inner and outer subcomponents. For comparison, the black histogram in each panel represents the system-to-system 
scatter in the Vrot/c for the stellar disc. The accreted, metal-poor, old, and/or outer spheroid subcomponents show only a mild degree of 
rotation and are usually dispersion-supported. The in situ, metal-rich, young, and/or inner spheroid subcomponents show higher rotation 
velocities on average but there is a large system-to-system scatter in Vrot/o" for these subcomponents. 



licity, age, or galactocentric distance. These are meant to 
roughly mimic what could be done observationally to dis- 
tinguish between the in situ and accreted components. On 
average, the stars that are metal poor, old, or preferentially 
located at large radii show only a mild degree of rotation 
and are usually dispersion supported, as was the case for 
the accreted component shown in the top left panel. Stars 
that are metal rich, young, or located in the inner regions 
of the spheroid show higher rotation velocities on average, 
in accordance with the in situ component shown in the top 
left panel. 



4-2.2 Rotation - Internal view 

To facilitate comparisons with kinematic studies of the 
spheroid of the Milky Way, which are generally limited to 
the Solar neighbourhood, we show in Fig.[TT]the Toomre di- 
agram for star particles situated within cylindrical galacto- 
centric radii of 7.0 kpc and 10.0 kpc (the Sun is at » 8.5 kpc 
in ou r Galaxy; e.g.. lVanhollebeke et al.|[2009l : [Gillessen et al.l 
I2OO9I ) and within 2 kpc of the disc midplane {\Z\ < 2 kpc). 
The velocity coordinates, V, U, and W correspond to the 
tangential, radial, and Z velocities in the rest frame of the 
Solar neighbourhood. The rest frame for each of the sim- 
ulated disc galaxies is computed by taking the medians of 
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Figure 11. Stacked Toomre diagram for star particles in the 'So- 
lar neighbourhood' of the simulated disc galaxies. The contours 
enclose 50% of the particles for any given subcomponent. Left of 
the vertical dotted line motion is retrograde with respect to the 
disc. Relatively metal-rich and young stars (which are predomi- 
nantly formed in situ) always rotate in prograde fashion and with 
a lag that ranges from 0.0 to -0.7 times the rotation velocity of 
the disc (see text). Relatively metal-poor and old stars (predom- 
inantly accreted) show a much larger range of velocities with a 
non-negligible fraction on retrograde orbits. 



each of the three velocity components for cold star particles 
(i.e., those assigned to the disc) within the selected volume. 
To account for the fact that our simulated galaxies span a 
range of disc (and total) masses, we have elected to nor- 
malise the velocities by the median rotation velocity of disc 
particles in the Solar neighbourhood, Vlsr- To derive the 
average result from the simulations, we stack the Toomre 
diagrams computed for each individual galaxy. To the left 
of the dotted vertical line motion is retrograde with respect 
to the disc. 

In the top panel of Fig. [11] we compare the velocity 
distributions for the in situ and accreted components of the 
spheroid in the Solar neighbourhood. The solid red and blue 
contours enclose 50% of the particles in the in situ and ac- 
creted components, respectively. The solid black contour en- 
closes 50% of the star particles in the disc component. Since 
the velocity rest frame is computed from disc star particles, 
the disc component is situated &t V — Vlsr ~ 0. The in 
situ component of the spheroid shows a clear prograde rota- 
tion signature with a lag ranging from -0.7 to 0.0 (median is 
« —0.35) times the rotation velocity of the disc. For a disc 
galaxy with a rotation velocity similar to M31, ~ 240 km/s, 
the median lag would correspond to « —84 km/s. As noted 
above this is similar to the velocity lag reported for the large 
disc-like structure around M31 by llbata et al.l |2005') . The 
accreted component of the spheroid, by contrast, shows a 
much broader range of velocities, with a non- negligible frac- 
tion of stars on retrograde orbits with respect to the disc. 

In the middle and bottom panels of Fig. [TT] we use 
metallicity and age, respectively, to effectively select the in 
situ and accreted components using observable quantities. 
As expected, young and/or metal- rich stars have a distribu- 
tion that is very similar to that of the pure in situ selection 
in the top panel. Old and/or metal-poor selection yields a 
distribution very similar to that of the pure accreted com- 
ponent selection. 

Based on the above, we find that there is a significantly 
different velocity distribution for the in situ (relatively young 
and metal-rich) and accreted (relatively old and metal-poor) 
components in the 'Solar neighbourhood' of the simulated 
disc galaxies. Encouragingly, there is mounting observa- 
tional evidence for multiple components in the kinemat- 
ics of stars in our o wn Solar neighbourhood. Fo r example, 
ICarollo et al] l|2010l ) (see also iBeers et al.l 120111 ) find kine- 
matic evidence for two components, one relatively metal-rich 
and one relatively metal-poor. (We disc uss these intriguing 
results in more detail below.) Likewise, iNissen fc Schusten 
(2010) (see also iNissen fc Schuster! 1201 ll ) have used a sam- 
ple of 94 dwarf stars with high-resolution spectral measure- 
ments and found evidence for two components in the kine- 
matics when they split their sample into stars with rela- 
tively high or low [a/Fe] for their [Fe/H]. They concluded 
that the "high-a" stars origin was likely in situ. Qualita- 
tiv ely our simulations are co nsistent with the trends found 
by INissen fc Schustej l|2010l ), although it appears that the 
overlap between the in situ and accreted components in 
Fig. [TT] is lar ger than that for the two components identi- 
fied bv .Nissen fc Schuster] (|2010l ) (see their Fig. 3). We point 
out, however, that Fig.[TT]is the average (stacked) result and 
that we expect there to be a large range of behaviours for 
individual galaxies. 
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Figure 12. Scatter plots of the best-fit siiape parameters (flattening c/a and power-law index a; see Fig. [S]!, rotation velocity (Vrot), 
and in situ mass fraction of the sphe roid (see text) . For co mparison we show recent observational estimates for the shape and rotation 
of the Milky Way spheroid made by iDeason et all ll2011a| [b[l. The top panel is the same as Fig. |8]but we include it for completeness. 
The shape parameters, c/a and a, are clearly correlated with the rotation velocity and there is a particularly strong trend between the 
rotation speed and the degree of flattening of the spheroid. The in situ mass fraction is correlated with both the flattening and the 
rotation velocity, in the sense that systems with a large in situ mass fraction are typically flatter and rotating faster, although there is a 
large degree of scatter in these trends. 



4.3 Shape and Rotation together 

The shape and kinematics of a system are obviously not 
independent quantities, although we have treated them as 
such above. We now examine the connection between the 
two. 

In Fig. [T2]we show scatter plots of the best-fit shape pa- 
rameters (flattening c/a and power-law index a; see Fig. [8]), 
rotation velocity (K-ot), and in situ mass fraction. Note that 
the latter two quantities have been computed using the same 
star particle selection as in the fitting of parametric model to 
the mass distribution (i.e., we exclude star particles within 
\Z\ < 4 kpc and beyond a 3D radius rso > 40 kpc). Note 
also that the top panel is the same as Fig.[8]but we include it 



for completeness. We comment on origin of the relationship 
between c/a and a below. 

The two panels in the middle row of Fig. [T^] show how 
c/a and a are related to the rotation of the spheroid. There 
is a strong correlation between the rotation velocity of the 
spheroid and its mass structure - systems with a high de- 
gree of rotation tend to be highly flattened and their mass 
density tends to drop off more slowly with radius than for 
systems with a lower degree of rotation. The correlation be- 
tween the rotation velocity and flattening is even stronger 
than the one between the rotation velocity and the power- 
law index. Like the structural parameters, there is a large 
range in rotation velocities across our sample, with some 
spheroids rotating as fast as 150 km/s while others are very 
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Figure 13. Scatter plot of the ratio of tiie rotation velocity to 
the 3D velocity dispersion of the spheroid against the ellipticity 
of the stellar spheroid. The solid curves correspond to the prcdic- 
tions of the tensor virial t heorem applied to oblate spheroids (see 
iBinnev fc TremaindfigSTl. eqn. 4-95) with anisotropy parameters 
5 (= 1 - o-L/o"D = 0-0 > 0-l> 0-2, 0.3, 0.4, and 0.5 (top to bot- 
tom). The inverted orange triangle represents the median of the 
simulated spheroids. The simulated galaxies have a wide range of 
behaviours: for some the shape is set entirely by rotation (5 = 0) 
while for others anisotropy plays as large a role as rotation. The 
median corresponds to (5 0.2. The Milky Way has e ft: 0.3 — 0.5 
and Vrot/o" J; 0.5 (see text), implying <5 0.2 — 0.4. 



mildly counter-rotating with respect to the disc. The median 
rotation velocity of the spheroid is -1-38 km/s. 

Is there a direct causal link between the rotation and 
the shape of the spheroid? It is possible, for example, that 
velocity anisotropy is the true underlying cause of the flat- 
tened shape rather than rotation (as is the case with many 
giant elliptical galaxies). We attempt to shed some light on 
this in Fig. 1131 which shows a scatter plot of the average 
ratio of the rotation velocity to the 3D velocity dispersion 
of the spheroid against its ellipticity, e = 1 — c/a. The solid 
curves correspond to the predictions of the te nsor virial the- 
orem applied to oblate s pheroids (see, e.g., iBinnevI [l973 : 
iBinnev fc Tremaind [l987l . see eqn. 4-95 of the latter) with 
anisotropy parameters 5 (= 1 — crL/f^a;) = 0.0 , 0.1, 0.2, 0.3, 
0.4, and 0.5 (top to bottom). Note that the predictions are 
independent of how rapidly the galaxy rotates and, if the 
isodensity surfaces are similar concentric ellipsoids, they are 
also independent of the radial density profile of the spheroid. 

It is readily apparent that the spheroids of the simu- 
lated galaxies display a fairly wide range of behaviours, from 
being completely rotation-dominated (i.e., points clustered 
around the top curve) to having a significant contribution 
to the flattening from velocity anisotropy with 5 ranging as 
high as ~ 0.5. The median implied value of 5 ~ 0.2 for the 
simulated galaxies. Thus, typically, both rotation and veloc- 
ity anisotropy contribute to the flattening of the spheroids 
in our simulations. 



Measuring the rotation speed of the stellar halo of the 
Milky Way is non-trivial as it relies on the availability of 
accurate distances to the spheroid tracers, an accurate dis- 
tance from the Sun to the centre of the Galaxy, and finally 
an accurate estimate of the local standard of rest circular 
velocity. While the Sun's galactocentric distance is now rel- 
atively well pinned down (to 8.5 kpc), the other two issues 
remain. ICarollo et all l|2007l . [2OI0I ) have used a large sample 
of main sequence turn off stars in the SDSS to constrain the 
rotational velocity of both the metal-rich inner and metal- 
poor outer components (both confined to rsD < 40 kpc). 
Under the assumption of the canonical local standard of rest 
circular velocity of 220 km/s, these authors find no evidence 
for rotation of the inner metal-rich component, but they find 
evidence for significant net retrograde motion of the outer 
metal-poor component of ~ —80 km/s. We note that there 
is an ongoing debate about whether or not hard kinematic 
evidence exists indicating the need for at least two halo com- 
po nents. ? argue t hat th e evidence for two components found 
bv lCaroUo etai] l|2010l ) is spurious, and the result of errors 
in the distan ce estimates t hat w ere employed in that study. 
However, as iBeers et al. I (I2OIII) demonstrate, the stars in 
the sample considered bv ICaroUo' et al. (201y) (whether or 
not the unlikely main-sequence turn-off luminosity classifi- 
cations are reassigned to dwarf and/or giant classifications) 
clearly exhibit proper motion distributions that are asym- 
metri c with respect to the other (non highly retrograde) 
stars. iBeers et al.1 l|201ll ) conclude this illustrates that the 
inferred asymmetry in the highly retrograde tail arises not 
from improper distances, but rather, is due primarily to real 
differences in the proper motions. 

iDeason et all l|2011ah have measured the rotation of the 
halo using a spectroscopic sample of BHB stars selected from 
the SDSS. The advantage of BHB stars is that they are lumi- 
nous and have nearly a constant absolute magnitude when 
selected within a certain colour range. The disadvantages 
are that there are far fewer BHB stars than main sequence 
turn off stars to work with, and that the BHBs may be a 
biased tracer of the spheroid (e.g., they are typically older 
on averag e than main se quence stars) . In qualitative agree - 
ment with lCaroUo et al.1 (2007, 2 01oh . lDeason et al.1 (|2011ah 
find evidence for a retrograde-rotating, relatively metal-poor 
component with Vrot ~ —25 km/s. But they also find evi- 
dence for a prograde-rotating, relatively metal-rich compo- 
nent with Vrot ~ 15 ± 8 km/s. These values were derived 
under the assumption of the canonical local standard of rest 
circular velocity of 22 km/s. A number o f recent studies 
(see the discussion in iDeason et al.ll2011al ) have, however, 
recommended an up ward revision to « 240 — 250 km/s. 
IDeason et al.l (|2011al ) show that if one adopts a local stan- 
dard of rest circular velocity of 240 km/s, then the impli- 
cation is that the metal-poor component has no detectable 
net rotation while the comparatively metal-rich component 
has a net prograde rotation velocity of « 45 km/s. This is 
the speed that is indicated by the solid blue cross in the 
two panels in the middle row of Fig. 1121 which is compatible 
with the distribution of simulated galaxies. 

Is the shape of the Milky Way's stellar halo determined 
by rotation or velocity anisotropy? We can use the theoret- 
ical curves plotted in Fig. [T3] to try to answer this ques- 
tion. As discussed in Section 4.1.2, the inferred shape of the 
Milky Way's halo is c/a ~ 0.5 — 0.7, implying an elliptic- 
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ity of e ~ 0.3 — 0.5. Using the latest BHB-based results of 
iDeason et al.l (l2011a ) provi des a measure of the rotation ve- 
locitv and lXue et al l|2008h have estimated the line-of-sight 
velocity dispersion of a similar BHB sample (see their Fig. 
10). But note that what is plotted in Fig. [13] is the ratio of 
the rotation velocity to the 3D velocity dispersion, which is 
not directly observable. Thus, one must assume something 
about the degree of anisotropy in order to infer an anisotropy 
from Fig. 1131 which is clearly not self-consistent. However, it 
can be said with some certainty from the BHB results that 
Kot/o" < 0.5 for the Milky Way's stellar halo. Using the 
theoretical curves in Fig. 1131 implies S « 0.2 — 0.4, which is 
consistent with the implied spread in the spheroids of the 
simulated galaxies. Thus, based on recent BHB measure- 
ments (which show the halo has a non-zero rotation speed), 
one infers that both rotation and velocity anisotropy con- 
tribute to the flattening of the Milky Way's stellar halo. 

We can now return to the top panel of Fig. 1121 to try to 
explain the relationship between c/a and a in terms of the 
correlations of these parameters with the rotational speed, 
Kot. There is a clear lower envelope to the distribution in 
the top panel: the minimum c/a increases with a (or, equiv- 
alently, the maximum a increases with c/a). This same en- 
velope is also visible in the middle left panel: the maximum 
a increases with decreasing Vrot and the central panel shows 
that c/a increases with decreasing Kot. Thus, the maximum 
a must increase with increasing c/a. However, this kind of 
correlation argument does not really explain the trend be- 
tween c/a and a or demonstrate causality. A possible expla- 
nation may be found in the Viot — ot plane. In particular, for a 
given a, Kot cannot be arbitrarily high. If it is too high, the 
centripetal force would exceed the gravitational one (which 
depends on a). 

Finally, in the bottom row of Fig. [T2] we investigate 
whether the shape of the stellar mass distribution or its 
kinematics tell us about the in situ mass fraction. The green 
points and error bars represent the median and 1-sigma (i.e., 
encloses 68% of the data) in the in situ mass fraction in bins 
of a, c/a, and Kot. As noted above, the in situ mass frac- 
tion has been computed within the volume limits specified 
above. 

There is a wide system-to-system variation in the in situ 
mass fraction, with typical fractions in the range of 20-60% 
(median of ~ 40%). We find no strong trend between the 
in situ mass fraction and the power-law index, a. A trend 
exists between the in situ mass fraction and the flattening 
of the spheroid for c/a < 0.6 or so, in that systems with 
large in situ mass fractions are typically flatter than those 
with low in situ mass fractions. An even stronger trend is 
evident between the in situ mass fraction and the rotation 
velocity of the stellar spheroid. If one uses the simulations 
as a guide to infer an in situ mass fraction for the Milky 
Way, the shape and kinematics would be consistent with a 
mass fraction of ~ 40%. 



5 SUMMARY 

We have used the GIMIC simulations to study the global 
structure and kine matics of the sphe roids of 412 simulated 
~ L* disc galaxies. iFont et al] l|201ll ) have recently demon- 
strated that these simulations are able to successfully repro- 



duce the metallicity and surface brightness proflles of the 
spheroids of the Milky Way and M31. A key to the success of 
the simulations is a significant contribution to the spheroid 
from stars that formed in situ. While the outer halo is domi- 
nated by accreted stars, stars formed in the main progenitor 
of the galaxy dominate at r < 30 kpc. 

In the present study we have shown that there is a close 
physical connection between the in situ spheroid and the 
disc, namely that the in situ spheroid was produced by the 
dynamical heating of a proto-disc. The likely mechanism 
for the heating is accretion of satellites (and mass in gen- 
eral) at z ~ 1 . A sim ilar conclusion was reached recently by 
iPurcell et af] (|2010l ) on the basis of the analysis of idealised 
simulations of disc bombardment by infalling satellites. 

Given the very different origins for the in situ and ac- 
creted spheroids, we expect there to be differences in the 
global structure and kinematics of the two components. In- 
deed, we have shown that the inner regions of the spheroid 
are highly flattened (oblate) and show evidence for two com- 
ponents in the kinematics, with the dominant one rotating 
in prograde fashion with respect to the disc. The degree 
of flattening is in good accordance with detailed observa- 
tions of both local galaxies, such as t he Milky Way (e.g 



Chiba fc Beer3l200ll : |jun'c et al.ll2008l:lDeason et al 



2011b 



Sesar et al, .1120111)" and M31 (jPritchet fc van den Berghlll994l : 



Ibata et al.ll2005 ), and stacked observations of more distant 
galaxies feib etti et al. liooi). This is in stark contrast to 
the results of collisionless hybrid models in which star parti- 
cles are 'paint ed on' to dark matter haloes of infalling dwarf 
galaxies (e.g., [B ullock fc Johnst"onll2005l: iFont et al.l [200631 : 
iDe Lucia fc He lmi 2008: ICooper et al.ll20ld '). which usually 
give rise to highly prolate c onfigurations, like t hat of the 
underlying dark matter halo jCooper et al.|[2010l ). 

As we noted in the Introduction, hybrid models also 
fail to reproduce the observed radial tr end with metallic- 
ity found in the Milky Way and M31 (iFont et all l2006bl : 
iDe Lucia fc Helmill2008l : ICooper et"alll2010l '). and it is dif- 
ficult to see how they could be easily reconciled with 
the kinematic and chemodyna mic evidence for a 'dua l 
halo' in the M ilky Way (e.g., ICarollo et all |2007| . l2010l : 
iNissen fc Schus ter 2010). The GiMiC simulations, on the 
other hand, successfully and simultaneously account for 
the large-scale metallicity distribution, the flattening of the 
spheroid, and the presence of two components in the kine- 
matics of the spheroid stars. 

Including hydrodynamics in the simulations does not 
in and of itself guarantee success though. We suggest that 
the stars that constitute the present-day in situ spheroid 
must be formed in a proto-disc (at relatively late times of 
z ~ 1 — 1.5) in order for the present-day spheroid to have the 
correct shape and kinematics. This likely requires an efficient 
form of feedback to operate in order to delay the formation 
of these stars. This is achieved in the GiMiC simulations with 
efficient, but energetically feasible, supernova feedback. Note 
that strong SN feedback is required in simulations in order 
to explain the generally low star formation efficiency and the 
'missing baryon problem' of normal galaxies, and to enrich 
the intergalactic medium via outflows. 

As a further test, it will be interesting to see whether 
such hydrodynamic simulations can also account for the rela- 
tively low fraction of unbound substructure (relative to that 
predicted by hybrid models) that is observed in the Milky 
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Way (e.g., IXue et al]|201lh . We intend to explore this in a 
future study. 
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APPENDIX: NUMERICAL CONVERGENCE 

Here we investigate the sensitivity of our results to numerical 
resolution. As noted in Section 2, the ~2a sphere has been 
simulated at eight times better mass resolution than the 
intermediate-resolution runs analysed in the present study. 

In Fig. [14] we compare the mean structure and kinemat- 
ics of 50 Milky Way- mass disc galaxies in the high-resolution 
simulation with those of the 412 galaxies in the intermediate- 
resolution runs. Overall, the agreement between the high- 
resolution and intermediate-resolution GIMIC simulations is 
quite good for both mean structure and kinematics. 
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Figure 14. Left: Stacked SDSS i-band surface brightness contour maps of the simulated galaxies, as in Fig. O for the intermediate- 
resolution (solid curves) and high-resolution (dashed curves) simulations. Right: Stacked line of sight {vy) velocity contour maps, as 
in Fig. |9] for the intermediate-resolution (solid curves) and high-resolution (dashed curves) simulations. Bluer colours indicate motion 
towards the observer and redder colours indicate motion away from the observer. Overall there is very good agreement between the mean 
structure and kinematics of the intermediate- and high-resolution galaxy samples. 
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